
%Take draws from the parameter space. 
%These draws are held in BETAS which is size NPxNVxNDRAWS
%Then calculate the probability of each person's sequence of choices at
%each draw for that person. 
%These probabilities are held in PROBS which is size NPxNDRAWS

BETAS=randi(NGridPts,[NP,NV,NDRAWS]); %BETAS go from 1 to NGridPts
BETAS=(BETAS-1)./(NGridPts-1);   %Now BETAS go from zero to 1 inclusive
for r=1:NV
  BETAS(:,r,:)=COEF(r,1)+(COEF(r,2)-COEF(r,1)).*BETAS(:,r,:);  %Now BETAS go from lower limit to upper limit for each coefficient
end;

rrate = squeeze(BETAS(:,1,:));
QH = squeeze(BETAS(:,2,:));
Psi = zeros(NP,T+1,NDRAWS);

Psi(:,1,:) = 1;
for ii=1:T
  Psi(:,(ii+1),:)= QH .* ((1 ./(1+rrate)).^(ii));
end
    
%expand to correct number of rows using the ID column
Psi = Psi(XMAT(:,1),:,:);

% load values of discounted variables
%FirstRun will change to zero in bootstrapping, so rows align with
%bootstrap sample
if (FirstRun==1 && subset==0)
    pricetime = csvread('price.time.csv', 1)./10;
    captime = csvread('CAP.time.csv', 1);
    swtime = csvread('SW.time.csv', 1);
    buffertime = csvread('buffer.time.csv', 1);
    jobstime = csvread('jobs.time.csv', 1);
    infratime = csvread('infra.time.csv', 1);
    wlifetime = csvread('wlife.time.csv', 1);
    qualitytime = csvread('quality.time.csv', 1);
elseif (FirstRun==1 && subset==1)
    pricetime = csvread('price.time1.csv', 1)./10;
    captime = csvread('CAP.time1.csv', 1);
    swtime = csvread('SW.time1.csv', 1);
    buffertime = csvread('buffer.time1.csv', 1);
    jobstime = csvread('jobs.time1.csv', 1);
    infratime = csvread('infra.time1.csv', 1);
    wlifetime = csvread('wlife.time1.csv', 1);
    qualitytime = csvread('quality.time1.csv', 1);
elseif (FirstRun==1 && subset==2)
    pricetime = csvread('price.time2.csv', 1)./10;
    captime = csvread('CAP.time2.csv', 1);
    swtime = csvread('SW.time2.csv', 1);
    buffertime = csvread('buffer.time2.csv', 1);
    jobstime = csvread('jobs.time2.csv', 1);
    infratime = csvread('infra.time2.csv', 1);
    wlifetime = csvread('wlife.time2.csv', 1);
    qualitytime = csvread('quality.time2.csv', 1);
elseif (FirstRun==1 && subset==3)
    pricetime = csvread('price.time3.csv', 1)./10;
    captime = csvread('CAP.time3.csv', 1);
    swtime = csvread('SW.time3.csv', 1);
    buffertime = csvread('buffer.time3.csv', 1);
    jobstime = csvread('jobs.time3.csv', 1);
    infratime = csvread('infra.time3.csv', 1);
    wlifetime = csvread('wlife.time3.csv', 1);
    qualitytime = csvread('quality.time3.csv', 1);    
elseif (FirstRun==1 && subset==4)
    pricetime = csvread('price.time4.csv', 1)./10;
    captime = csvread('CAP.time4.csv', 1);
    swtime = csvread('SW.time4.csv', 1);
    buffertime = csvread('buffer.time4.csv', 1);
    jobstime = csvread('jobs.time4.csv', 1);
    infratime = csvread('infra.time4.csv', 1);
    wlifetime = csvread('wlife.time4.csv', 1);
    qualitytime = csvread('quality.time4.csv', 1);  
end

dbenefits = zeros(NROWS, 7, NDRAWS);
dprice = zeros(NROWS, NDRAWS);
for ii=1:NDRAWS
        dprice(:,ii) = sum(squeeze(Psi(:,:,ii)) .* pricetime, 2);
        dbenefits(:,1,ii) = sum(squeeze(Psi(:,:,ii)) .* captime, 2);
        dbenefits(:,2,ii) = sum(squeeze(Psi(:,:,ii)) .* swtime, 2);
        dbenefits(:,3,ii) = sum(squeeze(Psi(:,:,ii)) .* buffertime, 2);
        dbenefits(:,5,ii) = sum(squeeze(Psi(:,:,ii)) .* jobstime, 2);
        dbenefits(:,4,ii) = sum(squeeze(Psi(:,:,ii)) .* infratime, 2);
        dbenefits(:,7,ii) = sum(squeeze(Psi(:,:,ii)) .* wlifetime, 2);
        dbenefits(:,6,ii) = sum(squeeze(Psi(:,:,ii)) .* qualitytime, 2);
end

%Add up the non-price variables times their WTP
v=sum(bsxfun(@times,dbenefits,BETAS(XMAT(:,1),4:end,:)),2);      %sum(NROWSx(NV-1)xNDRAWS,2) gives NROWSx1xNDRAWS
v=squeeze(v);                              %NROWSxNDRAWS
v=bsxfun(@minus,v,dprice); %Subtract price
v=squeeze(BETAS(XMAT(:,1),3,:)).*v; %Multiply everything by price/scale coef
v=exp(v);
v(isinf(v))= 1.8e307;
sparsematrix=bsxfun(@eq,sparse(1:NCS)',XMAT(:,2)'); %NCSxNROWS
denom=double(sparsematrix)*v;                           %NCSxNDRAWS
denom(denom==0)=0.0000000000000001;
p=v(XMAT(:,3)==1,:)./denom;                     %NCSxNDRAWS
p(p==0)=0.0000000000000001;
personid=XMAT(XMAT(:,3)==1,1);                  %NCSx1
sparsematrix=bsxfun(@eq,sparse(1:NP)',personid');    %NPxNCS
PROBS=double(sparsematrix)*log(p);                      %NPxNDRAWS
PROBS=exp(PROBS);                               %NPxNDRAWS
clear sparsematrix v denom p personid dprice dbenefits Psi rrate
if WantBoot==0
    clear qualitytime wlifetime infratime jobstime buffertime swtime captime
end



